x <- 1:400
y <- sin(x / 10) * exp(x * -0.01)
plot(x,y)
philly <- readOGR("Philly3")
## OGR data source with driver: ESRI Shapefile
## Source: "Philly3", layer: "Philly3"
## with 381 features
## It has 12 fields
names(philly)
## [1] "OBJECTID" "AREA" "PERIMETER" "TRACTID" "FIPSSTCO"
## [6] "TRT2000" "STFID" "N_HOMIC" "TOTALPO" "mdHHnc_"
## [11] "HOMIC_R" "PCT_COL"
plot(philly)
spplot(philly)
# RColorBrewer for colors
library(RColorBrewer)
display.brewer.all(type="seq")
pal <- brewer.pal(5, 'OrRd')
typeof(pal) #character
## [1] "character"
philly$HOMIC_R
## [1] 426.74253 22.70148 122.45898 24.52784 22.70148 29.73536
## [7] 122.45898 65.25711 122.45898 65.25711 122.45898 122.45898
## [13] 406.12309 29.73536 122.45898 65.25711 122.45898 19.49698
## [19] 42.58944 65.25711 595.23810 0.00000 22.41147 34.19973
## [25] 0.00000 118.97680 0.00000 0.00000 0.00000 0.00000
## [31] 126.18297 595.23810 141.37606 72.93946 22.41147 46.88233
## [37] 0.00000 328.46715 0.00000 34.19973 118.97680 0.00000
## [43] 190.80328 22.41147 0.00000 0.00000 0.00000 0.00000
## [49] 190.80328 66.22517 328.46715 66.22517 46.88233 59.48840
## [55] 19.49698 72.93946 190.80328 141.37606 595.23810 0.00000
## [61] 259.74026 141.37606 66.22517 0.00000 269.98650 259.74026
## [67] 0.00000 66.22517 66.22517 193.67334 259.74026 72.93946
## [73] 259.74026 48.94763 95.58402 259.74026 299.53917 595.23810
## [79] 66.22517 109.52903 40.76641 0.00000 259.74026 211.41649
## [85] 40.76641 86.13264 193.67334 88.82309 259.74026 95.58402
## [91] 595.23810 88.82309 0.00000 126.28541 86.13264 595.23810
## [97] 259.74026 180.73845 95.58402 180.73845 211.41649 88.82309
## [103] 95.58402 245.48092 109.52903 595.23810 72.96607 111.11111
## [109] 86.13264 88.82309 40.76641 279.13468 70.12623 372.12449
## [115] 175.39464 100.83549 180.73845 134.37850 109.52903 45.46143
## [121] 465.11628 121.40132 86.13264 82.06138 0.00000 294.77821
## [127] 279.13468 100.83549 217.83526 45.46143 180.73845 0.00000
## [133] 45.46143 877.54288 877.54288 217.83526 253.36500 279.13468
## [139] 279.13468 502.18694 502.18694 180.73845 0.00000 0.00000
## [145] 45.46143 330.50047 45.46143 0.00000 285.19575 0.00000
## [151] 260.75619 306.69620 638.29787 285.19575 330.50047 877.54288
## [157] 253.67834 638.29787 45.46143 918.03279 502.18694 0.00000
## [163] 350.20694 877.54288 226.78594 285.19575 260.75619 877.54288
## [169] 0.00000 1343.87352 285.19575 226.78594 350.20694 663.50711
## [175] 136.96285 638.29787 328.52450 404.70225 244.89796 638.29787
## [181] 23.18572 136.96285 350.20694 663.50711 226.78594 656.25707
## [187] 117.78563 988.21741 244.89796 244.89796 663.50711 656.25707
## [193] 226.78594 288.73917 112.88805 426.49272 663.50711 456.10034
## [199] 240.32684 57.58157 456.10034 155.35992 567.06114 240.32684
## [205] 567.06114 112.88805 663.50711 53.03633 112.88805 663.50711
## [211] 663.50711 537.41215 537.41215 988.21741 755.12406 689.51931
## [217] 293.82958 302.43729 726.16538 726.16538 38.66976 293.82958
## [223] 1080.19640 539.32584 1067.86083 223.71365 237.60494 210.18216
## [229] 38.66976 726.16538 539.32584 54.79452 539.32584 53.03633
## [235] 38.66976 538.07272 203.09723 781.44694 764.08787 293.82958
## [241] 449.16292 210.18216 781.44694 163.20821 479.23323 781.44694
## [247] 54.79452 465.11628 347.03684 516.94428 293.82958 553.86320
## [253] 764.08787 137.42556 479.23323 210.18216 28.54696 600.24010
## [259] 347.03684 465.11628 465.11628 764.08787 764.08787 449.16292
## [265] 421.11173 236.62177 449.16292 764.08787 348.04909 347.03684
## [271] 764.08787 553.86320 137.42556 68.44627 280.11204 449.16292
## [277] 68.44627 301.84123 280.11204 245.26980 553.86320 240.21963
## [283] 354.05193 347.03684 312.06657 613.49693 613.49693 301.84123
## [289] 220.88987 663.34992 663.34992 220.88987 312.06657 287.52156
## [295] 312.06657 470.58824 392.43667 159.13431 380.16791 470.58824
## [301] 177.17034 177.17034 250.11034 27.70851 0.00000 0.00000
## [307] 38.41721 38.41721 27.70851 143.64376 143.64376 0.00000
## [313] 143.64376 177.17034 0.00000 104.42047 143.64376 250.11034
## [319] 345.78147 147.71049 354.47026 147.71049 40.46945 175.47280
## [325] 104.42047 354.47026 104.42047 74.57122 40.46945 343.81139
## [331] 74.57122 345.78147 175.47280 354.47026 74.57122 382.40918
## [337] 563.44377 354.47026 74.57122 74.57122 594.87179 880.00000
## [343] 115.63367 115.63367 608.64273 363.14117 354.47026 880.00000
## [349] 40.01601 117.39845 363.14117 240.26910 354.47026 240.26910
## [355] 880.00000 78.35455 54.47612 880.00000 40.01601 76.96254
## [361] 240.26910 54.47612 240.26910 0.00000 160.10247 124.22360
## [367] 360.03600 54.47612 0.00000 54.47612 54.47612 0.00000
## [373] 0.00000 304.41400 0.00000 304.41400 0.00000 0.00000
## [379] 304.41400 567.06114 567.06114
spplot(philly,"HOMIC_R", col.regions = pal, cuts = 4) #applying colors and cuts
#===============================================================
# Denver datasets
denver1 <- readOGR("county_boundary_lines")
## OGR data source with driver: ESRI Shapefile
## Source: "county_boundary_lines", layer: "county_boundary_lines"
## with 936 features
## It has 13 fields
plot(denver1)
# spplot(denver1)
denver2 <- readOGR("county_boundary")
## OGR data source with driver: ESRI Shapefile
## Source: "county_boundary", layer: "county_boundary"
## with 12 features
## It has 1 fields
plot(denver2)
spplot(denver2)
# denver3 <- readOGR("tl_2017_us_county")
# plot(denver3)
# spplot(denver3)
denver4 <- readOGR("statisticalneighborhoods")
## OGR data source with driver: ESRI Shapefile
## Source: "statisticalneighborhoods", layer: "statistical_neighborhoods"
## with 78 features
## It has 4 fields
plot(denver4)
spplot(denver4)
names(denver4)#"NBHD_ID" "NBHD_NAME" "TYPOLOGY" "NOTES"
## [1] "NBHD_ID" "NBHD_NAME" "TYPOLOGY" "NOTES"
denver4$NBHD_NAME
## [1] Chaffee Park Sunnyside
## [3] Highland Globeville
## [5] Jefferson Park Sun Valley
## [7] Valverde Athmar Park
## [9] Windsor Northeast Park Hill
## [11] Elyria Swansea Wellshire
## [13] University Rosedale
## [15] Cheesman Park Hilltop
## [17] Montclair Hale
## [19] North Park Hill South Park Hill
## [21] University Park Platt Park
## [23] College View - South Platte Overland
## [25] Ruby Hill Kennedy
## [27] Hampden Baker
## [29] Fort Logan Bear Valley
## [31] Harvey Park South Southmoor Park
## [33] Hampden South Indian Creek
## [35] Goldsmith Virginia Village
## [37] Gateway - Green Valley Ranch DIA
## [39] University Hills Harvey Park
## [41] Mar Lee Westwood
## [43] East Colfax Auraria
## [45] Cory - Merrill Belcaro
## [47] Washington Park Washington Park West
## [49] Speer Cherry Creek
## [51] Country Club Congress Park
## [53] City Park Clayton
## [55] Skyland Cole
## [57] Marston Washington Virginia Vale
## [59] Barnum Barnum West
## [61] Villa Park West Colfax
## [63] West Highland Sloan Lake
## [65] Berkeley Regis
## [67] Lincoln Park City Park West
## [69] Whittier Capitol Hill
## [71] North Capitol Hill Civic Center
## [73] CBD Union Station
## [75] Five Points Stapleton
## [77] Montbello Lowry Field
## 78 Levels: Athmar Park Auraria Baker Barnum Barnum West ... Windsor
# denver4$TYPOLOGY #NA
spplot(denver4, "NBHD_NAME", col.regions = pal, cuts = 78)
d5pal <- brewer.pal(8, 'YlOrRd')
denver5 <- readOGR("animal_care_and_control_division_boundaries")
## OGR data source with driver: ESRI Shapefile
## Source: "animal_care_and_control_division_boundaries", layer: "animal_care_and_control_division_boundaries"
## with 8 features
## It has 4 fields
plot(denver5)
# spplot(denver5)
names(denver5)#"DISTRICT_I" "INSPECTOR_" "INSPECTOR1" "INSPECTO_1"
## [1] "DISTRICT_I" "INSPECTOR_" "INSPECTOR1" "INSPECTO_1"
denver5$DISTRICT_I
## [1] F26 F27 F29 F32 F31 F33 F30 F28
## Levels: F26 F27 F28 F29 F30 F31 F32 F33
# denver5$INSPECTO_1 #INSPECTO_1
spplot(denver5, "DISTRICT_I", col.regions = d5pal, cuts = 78)
library(classInt)
breaks_qt <- classIntervals(philly$HOMIC_R, n = 5) # 5 colors
str(breaks_qt)
## List of 2
## $ var : num [1:381] 426.7 22.7 122.5 24.5 22.7 ...
## $ brks: num [1:6] 0 53 134 261 471 ...
## - attr(*, "style")= chr "quantile"
## - attr(*, "nobs")= int 156
## - attr(*, "call")= language classIntervals(var = philly$HOMIC_R, n = 5)
## - attr(*, "intervalClosure")= chr "left"
## - attr(*, "class")= chr "classIntervals"
breaks_qt$brks
## [1] 0.00000 53.03633 134.37850 260.75619 470.58824 1343.87352
spplot(philly, "HOMIC_R", col.regions=pal, at = breaks_qt$brks, main = "Philadelphia homicide rate per 100,000")
# add a very small value to the top breakpoint, and subtract from the bottom for symmetry
br <- breaks_qt$brks
offs <- 0.0000001
br[1] <- br[1] - offs
br[length(br)] <- br[length(br)] + offs
# plot
spplot(philly, "HOMIC_R", col.regions=pal, at=br, main = "Philadelphia homicide rate per 100,000")
# correcting legend
philly$HOMIC_R_bracket <- cut(philly$HOMIC_R, br)
head(philly$HOMIC_R_bracket)
## [1] (261,471] (-1e-07,53] (53,134] (-1e-07,53] (-1e-07,53] (-1e-07,53]
## Levels: (-1e-07,53] (53,134] (134,261] (261,471] (471,1.34e+03]
class(philly$HOMIC_R_bracket)
## [1] "factor"
spplot(philly, "HOMIC_R_bracket", col.regions=pal, main = "Philadelphia homicide rate per 100,000")
# install.packages('GISTools')
library(GISTools) # load library
## Loading required package: MASS
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
# View(philly)
choropleth(philly, philly$HOMIC_R) # plot the polygons
shd <- auto.shading(philly$HOMIC_R) # we need that for the legend coloring
choro.legend( # plot the legend
bbox(philly)["x","max"] - 5000, # x coordinate of top left corner
bbox(philly)["y","min"] + 15000, # y coordinate of top left corner
shd # color scheme
)
title("Philadelphia homicide rate per 100,000") # plot the title.
library(sf)
## Linking to GEOS 3.6.1, GDAL 2.1.3, proj.4 4.9.3
# getwd
philly_sf <- st_read("/Users/josephhaymaker/Desktop/STAT571_Denver_Animal_Protection/Philly3/Philly3.shp")
## Reading layer `Philly3' from data source `/Users/josephhaymaker/Desktop/STAT571_Denver_Animal_Protection/Philly3/Philly3.shp' using driver `ESRI Shapefile'
## Simple feature collection with 381 features and 12 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 1739498 ymin: 457288.9 xmax: 1764018 ymax: 490539.6
## epsg (SRID): NA
## proj4string: +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs
plot(philly_sf)
## Warning: plotting the first 10 out of 12 attributes; use max.plot = 12 to
## plot all
## Warning in classInt::classIntervals(values, nbreaks, breaks): var has
## missing values, omitted in finding classes
plot(philly_sf$HOMIC_R) # this is a numeric vector
plot(philly_sf["HOMIC_R"])
hr_cuts <- cut(philly_sf$HOMIC_R, br)
plot(philly_sf["HOMIC_R"], main = "Philadelphia homicide rate per 100,000", col = pal[as.numeric(hr_cuts)])
legend(1760000, 471000, legend = paste("<", round(br[-1])), fill = pal)
# GGPLOT
#So if we wanted to map polygons, like census tract boundaries, we would use longitude and latitude of their vertices as our x and y values and geom_polygon() as our geometry.
library(broom)
philly_df <- tidy(philly)
## Regions defined for each Polygons
dim(philly_df) #31410 7
## [1] 31410 7
# head(philly_df, 30)
# View(philly_df)
# head(philly_df)
str(philly, max.level = 2)
## Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 381 obs. of 13 variables:
## ..@ polygons :List of 381
## ..@ plotOrder : int [1:381] 378 377 11 182 25 16 372 3 94 360 ...
## ..@ bbox : num [1:2, 1:2] 1739498 457289 1764018 490540
## .. ..- attr(*, "dimnames")=List of 2
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
philly$polyID <- sapply(slot(philly, "polygons"), function(x) slot(x, "ID"))
philly_df <- merge(philly_df, philly, by.x = "id", by.y="polyID")
head(philly_df)
## id long lat order hole piece group OBJECTID AREA PERIMETER
## 1 0 1758028 490539.6 1 FALSE 1 0.1 1 37863520 32091.95
## 2 0 1758067 490520.2 2 FALSE 1 0.1 1 37863520 32091.95
## 3 0 1758134 490486.9 3 FALSE 1 0.1 1 37863520 32091.95
## 4 0 1758294 490407.7 4 FALSE 1 0.1 1 37863520 32091.95
## 5 0 1758369 490370.2 5 FALSE 1 0.1 1 37863520 32091.95
## 6 0 1758451 490329.6 6 FALSE 1 0.1 1 37863520 32091.95
## TRACTID FIPSSTCO TRT2000 STFID N_HOMIC TOTALPO mdHHnc_ HOMIC_R
## 1 365 42101 036500 42101036500 3 703 NA 426.7425
## 2 365 42101 036500 42101036500 3 703 NA 426.7425
## 3 365 42101 036500 42101036500 3 703 NA 426.7425
## 4 365 42101 036500 42101036500 3 703 NA 426.7425
## 5 365 42101 036500 42101036500 3 703 NA 426.7425
## 6 365 42101 036500 42101036500 3 703 NA 426.7425
## PCT_COL HOMIC_R_bracket
## 1 1.78117 (261,471]
## 2 1.78117 (261,471]
## 3 1.78117 (261,471]
## 4 1.78117 (261,471]
## 5 1.78117 (261,471]
## 6 1.78117 (261,471]
View(philly_df)
library(ggplot2)
ggplot() + # initialize ggplot object
geom_polygon( # make a polygon
data = philly_df, # data frame
aes(x = long, y = lat, group = group, # coordinates, and group them by polygons
fill = cut_number(HOMIC_R, 5))) + # variable to use for filling
scale_fill_brewer("Homicide Rate", palette = "OrRd") + # fill with brewer colors
ggtitle("Philadelphia homicide rate per 100,000") + # add title
theme(line = element_blank(), # remove axis lines ..
axis.text=element_blank(), # .. tickmarks..
axis.title=element_blank(), # .. axis labels..
panel.background = element_blank()) + # .. background gridlines
coord_equal() # both axes the same scale
# ggplot(philly_sf) + geom_sf(aes(fill=HOMIC_R))
library(tmap)
names(philly)
## [1] "OBJECTID" "AREA" "PERIMETER"
## [4] "TRACTID" "FIPSSTCO" "TRT2000"
## [7] "STFID" "N_HOMIC" "TOTALPO"
## [10] "mdHHnc_" "HOMIC_R" "PCT_COL"
## [13] "HOMIC_R_bracket" "polyID"
tm_shape(philly) +
tm_polygons("HOMIC_R", style="quantile", title="Philadelphia \nhomicide rate \nper 100,000") +
tm_shape(philly) + tm_polygons("TOTALPO", style="quantile") +
tmap_mode("view")
## tmap mode set to interactive viewing
last_map()
## tmap mode set to interactive viewing
vignette("tmap-nutshell")
## starting httpd help server ...
## done
# ============================================
tm_shape(denver4) +
tm_polygons("NBHD_NAME", style="quantile", title="Denver county statistical neighborhoods")
# tmap_mode("view")
# last_map()
tm_shape(denver5) +
tm_polygons("DISTRICT_I", style="quantile", title="Philadelphia \nhomicide rate \nper 100,000")
New tutorial - link
library(rgeos)
library(maptools)
library(ggplot2)